# Carga de librerías
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.2 ✓ purrr 0.3.4
## ✓ tibble 3.0.3 ✓ dplyr 1.0.0
## ✓ tidyr 1.1.0 ✓ stringr 1.4.0
## ✓ readr 1.3.1 ✓ forcats 0.5.0
## ── Conflicts ────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
library(tidymodels)
## ── Attaching packages ────────────────────────────────────────── tidymodels 0.1.1 ──
## ✓ broom 0.7.0 ✓ recipes 0.1.13
## ✓ dials 0.0.9 ✓ rsample 0.0.8
## ✓ infer 0.5.3 ✓ tune 0.1.1
## ✓ modeldata 0.1.0 ✓ workflows 0.2.1
## ✓ parsnip 0.1.4 ✓ yardstick 0.0.7
## ── Conflicts ───────────────────────────────────────────── tidymodels_conflicts() ──
## x scales::discard() masks purrr::discard()
## x dplyr::filter() masks stats::filter()
## x recipes::fixed() masks stringr::fixed()
## x dplyr::lag() masks stats::lag()
## x yardstick::spec() masks readr::spec()
## x recipes::step() masks stats::step()
library(epiDisplay)
## Loading required package: foreign
## Loading required package: survival
## Loading required package: MASS
##
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select
## Loading required package: nnet
##
## Attaching package: 'epiDisplay'
## The following object is masked from 'package:yardstick':
##
## kap
## The following object is masked from 'package:scales':
##
## alpha
## The following object is masked from 'package:ggplot2':
##
## alpha
library(GGally)
## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2
library(corrr)
library(knitr)
library(ggrepel)
library(kableExtra)
##
## Attaching package: 'kableExtra'
## The following object is masked from 'package:dplyr':
##
## group_rows
library(here)
## here() starts at /home/daro/cdatos/TPs/EEA/TP2
El dataset de precios de inmuebles proviene de Properati. El mismo ya fue filtrado por los docentes y a su vez se encuentra particionado en subconjuntos de training y test.
train_ds <- read_csv(here("/ds/ar_properties_train.csv"))
## Parsed with column specification:
## cols(
## id = col_character(),
## l3 = col_character(),
## rooms = col_double(),
## bathrooms = col_double(),
## surface_total = col_double(),
## surface_covered = col_double(),
## price = col_double(),
## property_type = col_character()
## )
head(train_ds)
dim_desc(train_ds)
## [1] "[32,132 x 8]"
summary(train_ds)
## id l3 rooms bathrooms
## Length:32132 Length:32132 Min. :1.000 Min. :1.000
## Class :character Class :character 1st Qu.:2.000 1st Qu.:1.000
## Mode :character Mode :character Median :3.000 Median :1.000
## Mean :2.722 Mean :1.427
## 3rd Qu.:3.000 3rd Qu.:2.000
## Max. :8.000 Max. :5.000
## surface_total surface_covered price property_type
## Min. : 28 Min. : 26.00 Min. : 69500 Length:32132
## 1st Qu.: 46 1st Qu.: 42.00 1st Qu.:120000 Class :character
## Median : 65 Median : 58.00 Median :169000 Mode :character
## Mean : 80 Mean : 70.61 Mean :214991
## 3rd Qu.: 99 3rd Qu.: 85.00 3rd Qu.:260000
## Max. :320 Max. :265.00 Max. :950000
unique(train_ds$l3)
## [1] "Almagro" "Palermo" "Caballito"
## [4] "Villa Crespo" "Villa Urquiza" "Barracas"
## [7] "Nuñez" "Belgrano" "Floresta"
## [10] "Recoleta" "Colegiales" "Parque Chas"
## [13] "Barrio Norte" "Villa Devoto" "Villa Ortuzar"
## [16] "Velez Sarsfield" "Parque Chacabuco" "Villa Pueyrredón"
## [19] "Villa Real" "Once" "Flores"
## [22] "San Telmo" "Las Cañitas" "Villa Santa Rita"
## [25] "Villa del Parque" "Retiro" "Parque Centenario"
## [28] "Chacarita" "Abasto" "San Cristobal"
## [31] "Liniers" "Balvanera" "Boedo"
## [34] "Villa General Mitre" "Saavedra" "Parque Patricios"
## [37] "Boca" "Paternal" "Monserrat"
## [40] "San Nicolás" "Parque Avellaneda" "Centro / Microcentro"
## [43] "Puerto Madero" "Congreso" "Monte Castro"
## [46] "Mataderos" "Coghlan" "Versalles"
## [49] "Villa Lugano" "Catalinas" "Villa Luro"
## [52] "Constitución" "Pompeya" "Tribunales"
## [55] "Agronomía" "Villa Soldati" "Villa Riachuelo"
unique(train_ds$property_type)
## [1] "Departamento" "PH" "Casa"
apply(train_ds, 2, function(x) any(is.na(x)))
## id l3 rooms bathrooms surface_total
## FALSE FALSE FALSE FALSE FALSE
## surface_covered price property_type
## FALSE FALSE FALSE
tab1(train_ds$property_type, sort.group = "decreasing", cum.percent = TRUE)
## train_ds$property_type :
## Frequency Percent Cum. percent
## Departamento 28203 87.8 87.8
## PH 3120 9.7 97.5
## Casa 809 2.5 100.0
## Total 32132 100.0 100.0
tab1(train_ds$l3, sort.group = "decreasing", cum.percent = TRUE)
## train_ds$l3 :
## Frequency Percent Cum. percent
## Palermo 4890 15.2 15.2
## Belgrano 2736 8.5 23.7
## Almagro 2486 7.7 31.5
## Caballito 2393 7.4 38.9
## Villa Crespo 2095 6.5 45.4
## Recoleta 1824 5.7 51.1
## Barrio Norte 1359 4.2 55.3
## Villa Urquiza 1319 4.1 59.4
## Flores 916 2.9 62.3
## Nuñez 887 2.8 65.1
## Balvanera 883 2.7 67.8
## Villa Devoto 527 1.6 69.4
## Colegiales 523 1.6 71.1
## Villa del Parque 500 1.6 72.6
## San Telmo 488 1.5 74.2
## Saavedra 433 1.3 75.5
## San Cristobal 417 1.3 76.8
## Parque Centenario 404 1.3 78.1
## Floresta 366 1.1 79.2
## Puerto Madero 357 1.1 80.3
## Boedo 336 1.0 81.3
## Paternal 327 1.0 82.4
## San Nicolás 325 1.0 83.4
## Monserrat 310 1.0 84.3
## Chacarita 292 0.9 85.3
## Retiro 291 0.9 86.2
## Villa Pueyrredón 285 0.9 87.0
## Parque Chacabuco 282 0.9 87.9
## Barracas 280 0.9 88.8
## Villa Luro 264 0.8 89.6
## Liniers 260 0.8 90.4
## Mataderos 259 0.8 91.2
## Once 255 0.8 92.0
## Congreso 246 0.8 92.8
## Coghlan 240 0.7 93.5
## Las Cañitas 187 0.6 94.1
## Abasto 179 0.6 94.7
## Parque Patricios 171 0.5 95.2
## Monte Castro 168 0.5 95.7
## Constitución 139 0.4 96.2
## Villa Santa Rita 128 0.4 96.6
## Villa Lugano 113 0.4 96.9
## Versalles 113 0.4 97.3
## Villa Ortuzar 109 0.3 97.6
## Centro / Microcentro 105 0.3 97.9
## Boca 104 0.3 98.3
## Villa General Mitre 99 0.3 98.6
## Parque Chas 78 0.2 98.8
## Parque Avellaneda 78 0.2 99.0
## Pompeya 65 0.2 99.2
## Velez Sarsfield 56 0.2 99.4
## Villa Real 55 0.2 99.6
## Tribunales 54 0.2 99.8
## Agronomía 51 0.2 99.9
## Villa Riachuelo 13 0.0 100.0
## Villa Soldati 9 0.0 100.0
## Catalinas 3 0.0 100.0
## Total 32132 100.0 100.0
lm_all <- lm(price ~ rooms + bathrooms + surface_total + surface_covered + property_type + l3, data = train_ds)
tidy_lm_all <- tidy(lm_all, conf.int = TRUE)
tidy_lm_all
(intercept) (-109143.7839): es una inferencia teórica del modelo, no representa la realidad, sería una propiedad sin superficie, habitaciones, barrio, etc.
rooms (-4359.7289): Al ser negativo implica una quita promedio de ese valor en el precio de la propiedad al aumentar en 1 la cantidad de ambientes.
bathrooms (34367.2123): Indica en cuanto aumenta en promedio el valor de la propiedad por agregar un baño.
surface_total (890.5585): Indica en cuanto aumenta en promedio el valor de la propiedad por aumentas en 1 m2 la superficie total.
surface_covered (1497.9310): Indica en cuanto aumenta en promedio el valor de la propiedad por aumentas en 1 m2 la superficie cubierta.
property_typeDepartamento (91485.9410): El valor de β (91485.9410) indica cuanto aumenta la función de respuesta para un departamento en comparación de una casa.
property_typePH (46220.0449): El valor de β (46220.0449) indica cuanto aumenta la función de respuesta para un PH en comparación de una casa.
l3{Barrio} : Son 56 coeficientes que indican en cuanto aumenta o se disminuye en promedio el precio de una propiedad dependiendo del barrio en donde esta ubicada.
Observación: para cada coeficiente se asume que se mantienen constantes el restos de las covariables.
Analizando los p-valores de las variables dummies asociadas a property_type estas resultan estadísticamente significativas. Por otro lado, las variables dummies asociados a la variable l3 se reparten entre algunas significativas y otras no.
tidy(anova(lm_all))
Analizando la significatividad global de las variables l3 y property_type y observando que tienen un p-valor menor a 0.5 se puede concluir que ambas son signigicativas.
glance(lm_all)
ggplot(tidy_lm_all,
aes(estimate, term, xmin = conf.low, xmax = conf.high, height = 1)) +
geom_point(color = "blue", size=3) +
geom_vline(xintercept = 0, lty = 2, color = "black") +
geom_errorbarh(color = "red", size=1) +
theme_bw() +
labs(y = "Coeficientes β", x = "Estimación")
El r.squared en el resúmen global representa el \(Rˆ2\) o coeficiente de determinación que refleja la bondad del ajuste del modelo a la variable precio. Al ser un valor cercano a 1 nos está diciendo que el modelo lm_all ajusta muy bien dicha variable.
Por otro lado, el gráfico nos permite visualizar con mayor claridad los intervalos de confiaza de de los coeficientes de las variables dummies de l3. Muchas de ellas incluyen el cero, ergo estas variables podrían considerarse no significativas.
lm_wo_l3 <- lm(price ~ rooms + bathrooms + surface_total + surface_covered + property_type, data = train_ds)
tidy_lm_wo_l3 <- tidy(lm_wo_l3, conf.int = TRUE)
tidy_lm_wo_l3
tidy(
anova(lm_wo_l3)
)
Al quitar la variable l3 la cantidad de coeficientes se reduce considerablemente, no necesitamos graficar para indentificar que ningúno de los intervalos de confianza incluye el 0. Por otro lado todos los p-values son menores a 0.5 ergo las variables se pueden considerar significativas, esto tambien se ve confirmado por el resultado del Test F.
models = list(lm_all = lm_all, lm_wo_l3 = lm_wo_l3)
purrr::map_df(models, broom::glance, .id = "model")
El \(Rˆ2\) de lm_all es superior al del modelo lm_wo_l3 por lo tanto el modelo que contiene todas las covariables es el que mejor explica la variable la variabilidad del precio.
En el análisis anterior se observó que l3 no es un buen agrupador, algunas de sus dummies eran significativas pero muchas otras no. Vamos a proceder a crear una nueva variable neighborhood que agrupe a las propiedad teniendo en cuenta el precio del metro cuadrado.
Es necesario tambien agregar otra variable auxiliar que contenta el precio por metro para luego poder hacer el agrupamiento necesario:
La valores posibles de neighborhood van a ser:
train_ds = train_ds %>% mutate(
price_m_sq = price / surface_total,
neighborhood = factor(
case_when(
price_m_sq <= quantile(price_m_sq)[2] ~ "low_price",
price_m_sq > quantile(price_m_sq)[2] & price_m_sq <= quantile(price_m_sq)[4] ~ "medium_price",
price_m_sq > quantile(price_m_sq)[4] ~ "high_price"
)
)
)
train_ds
lm_neighborhoods <- lm(price ~ rooms + bathrooms + surface_total + surface_covered + property_type + neighborhood, data = train_ds)
tify_lm_neighborhoods <- tidy(lm_neighborhoods, conf.int = TRUE)
tify_lm_neighborhoods
Observaciones:
models = list( lm_all = lm_all, lm_wo_l3 = lm_wo_l3, lm_neighborhoods = lm_neighborhoods )
purrr::map_df(models, broom::glance, .id = "model") %>% arrange(adj.r.squared)
El modelo lm_neighborhoods supera al lm_all en el \(Rˆ2\) y por lo tanto tiene mayor explicabilidad sobre la variabilidad del precio de las propiedades.
Las variables surface_total y surface_covered están muy correlacionadas por lo que vamos a generar una nueva variable llamada surface_uncovered que representa la diferencia entre la superficie total y la cubierta.
train_ds = train_ds %>% mutate(surface_uncovered = surface_total - surface_covered)
train_ds
lm_surface_uncovered <- lm( price ~ rooms + bathrooms + surface_uncovered + surface_covered + property_type + price_m_sq + neighborhood, data = train_ds)
tidy_lm_surface_uncovered <- tidy(lm_surface_uncovered, conf.int = TRUE)
tidy_lm_surface_uncovered
Observaciones:
models = list( lm_all = lm_all, lm_wo_l3 = lm_wo_l3, lm_neighborhoods = lm_neighborhoods, lm_surface_uncovered = lm_surface_uncovered )
purrr::map_df(models, broom::glance, .id = "model") %>% arrange(adj.r.squared)
El modelo lm_surface_uncovered supera al lm_neighborhoods en el \(Rˆ2\) y por lo tanto tiene mayor explicabilidad sobre la variabilidad del precio de las propiedades.
plot(lm_surface_uncovered)
Del plot de Residual vs Fitted se puede observar q no se respeta la linealidad parece que no se respeta y por otro lado mientras aumentan las predicciones aumenta la heterocedasticidad, además existen varios puntos con un leverage alto que se pueden observar en el plot Residual vs leverage, donde al menos dos pueden ser outliers. Por último en el plot **Normal Q-Q* se observa que los extremos no se ajustan a la distribución teórica.
En conclusión, el modelo lm_surface_uncovered no comple con los supuestos del modelo lineal.
\[ log(price) = \beta_0 + \beta_1log(rooms) + \beta_2log(bathrooms) + \beta_3log(surface\_covered) + \beta_4property\_type + \beta_5neighborhood + \beta_6surface\_uncovered \]
lm_log = lm( log(price) ~ log(rooms) + log(bathrooms) + log(surface_covered) + property_type + neighborhood + surface_uncovered, data = train_ds )
tidy_ml_log <- tidy(lm_log, conf.int = TRUE)
tidy_ml_log
tidy(anova(lm_log))
Observaciones:
models = list( lm_log = lm_log, lm_surface_uncovered = lm_surface_uncovered )
purrr::map_df(models, broom::glance, .id = "model")
El modelo lm_surface_uncovered supera al lm_logs en el \(Rˆ2\) y por lo tanto tiene mayor explicabilidad sobre la variabilidad del precio de las propiedades.
plot(lm_log)
Del plot de Residual vs Fitted se puede observar que se repeta cierta linealidad y, por otro lado, mientras aumentan las predicciones aumenta la heterocedasticidad, además se remarcan algunos puntos en el plot: 3636, 7346, 3709. En Residual vs leverage vemos varios puntos con un leverage alto y se remarcan algunos puntos tambien. Por último em Normal Q-Q se observa que los extremos no se ajustan a la distribución teórica.
En conclusión, el modelo lm_log no comple con los supuestos del modelo lineal.
Vamos a crear dos nuevos modelos
\[ log(price) = \beta_0 + \beta_1log(surface\_total) + \beta_2log(bathrooms) + \beta_3log(surface\_covered) + \beta_4property\_type + \beta_5neighborhood + \beta_6surface\_uncovered \]
lm_log_v2 = lm( log(price) ~ log(surface_total) + log(bathrooms) + log(surface_covered) + property_type + neighborhood + surface_uncovered, data = train_ds )
tidy_lm_log_v2 <- tidy(lm_log_v2, conf.int = TRUE)
tidy_lm_log_v2
Observaciones:
plot(lm_log_v2)
Podemos ver como el modelo se comporta de forma muy similar al lm_log.
\[ log(price) = \beta_0 + \beta_1log(surface\_total) + \beta_2log(bathrooms) + \beta_3log(surface\_covered) + \beta_4neighborhood + \beta_5surface\_uncovered \]
lm_log_v3 = lm( log(price) ~ log(surface_total) + log(bathrooms) + log(surface_covered) + neighborhood + surface_uncovered, data = train_ds )
tidy_lm_log_v3 <- tidy(lm_log_v3, conf.int = TRUE)
tidy_lm_log_v3
Observaciones:
plot(lm_log_v3)
Podemos ver como el modelo se comporta de forma muy similar al lm_log.
Se seleccionan para evaluación los modelos logaritmicos ml_log, ml_log_v2 y ml_log_v3 y el modelo lm_surface_uncovered.
models = list(lm_surface_uncovered = lm_surface_uncovered, lm_log = lm_log, lm_log_v2 = lm_log_v2, lm_log_v3 = lm_log_v3 )
purrr::map_df(models, broom::glance, .id = "model") %>% arrange(adj.r.squared)
El modelo lm_log_v2 supera al resto de los modelos en el \(Rˆ2\) y por lo tanto tiene mayor explicabilidad sobre la variabilidad del precio de las propiedades. De todas formas hay que destacar que todos los modelos tienen un \(Rˆ2\) mayor a 0.91.
test_ds <- read_csv(here("/ds/ar_properties_test.csv"))
## Parsed with column specification:
## cols(
## id = col_character(),
## l3 = col_character(),
## rooms = col_double(),
## bathrooms = col_double(),
## surface_total = col_double(),
## surface_covered = col_double(),
## price = col_double(),
## property_type = col_character()
## )
test_ds = test_ds %>%
mutate(
price_m_sq = price / surface_total,
neighborhood = factor(
case_when(
price_m_sq <= quantile(price_m_sq)[2] ~ "low_price",
price_m_sq > quantile(price_m_sq)[2] & price_m_sq <= quantile(price_m_sq)[4] ~ "medium_price",
price_m_sq > quantile(price_m_sq)[4] ~ "high_price"
)
),
surface_uncovered = surface_total - surface_covered
)
test_ds
models_log = list(lm_log = lm_log, lm_log_v2 = lm_log_v2, lm_log_v3 = lm_log_v3 )
pred_train_log = map(.x = models_log, .f = augment)
map_dfr(
.x = pred_train_log,
.f = rmse,
truth = exp(`log(price)`),
estimate = exp(.fitted),
.id="modelo"
) %>% arrange(.estimate)
pred_train_surface_uncovered = augment(lm_surface_uncovered)
rmse(
data = pred_train_surface_uncovered,
truth = price,
estimate = .fitted
)
De de los 4 modelos evaluados en training se destaca lm_surface_uncovered por tener el menor valor de RMSE.
pred_test_log = map(
.x = models_log,
.f = augment,
newdata = test_ds
)
map_dfr(
.x = pred_test_log,
.f = rmse,
truth = price,
estimate = exp(.fitted),
.id="modelo"
) %>%
arrange(.estimate)
pred_test_surface_uncovered = augment(
lm_surface_uncovered,
newdata=test_ds
)
pred_test_surface_uncovered
rmse(
data = pred_test_surface_uncovered,
truth = price,
estimate = .fitted
)
De de los 4 modelos evaluados en testing se destaca lm_surface_uncovered por tener el menor valor de RMSE al igual que con el dataset de training.
De los 4 modelos evaluados el que se destaca por sobre el resto es lm_surface_uncovered que tiene un \(Rˆ2\) alto de 0.9211874 que asegura una alta explicabilidad del precio de las propiedades y además tiene el RMSE (error cuadrático medio) mas bajo lo cuál nos asegura que sus predicciones van a ser las mas cercanas a los valores reales. Es para destacar que los 4 modelos seleccionados tienen un \(Rˆ2\) mayor a 0.91 y que ninguno parecen cumplir los supuestos del modelo lineal.